library(spatstat)
library(tidyverse)

data<-read.csv("mcd_v_mgn_working.csv", header = TRUE, stringsAsFactors = FALSE,
               check.names = FALSE)

dup_data<-duplicated(data)
data<-data[!dup_data,]
rm(dup_data)

images<-unique(data$image)

###################################Windows######################################
tiles<- read.csv("tiles.csv", header = TRUE, stringsAsFactors = FALSE, 
                 check.names = FALSE)

mask_gen<-function(a, b, d, e){
  mask_list<- vector("list", length = length(a))
  for (i in 1:length(a)){
     mask_df<-as.data.frame(matrix(ncol = 2))
     colnames(mask_df)<- c("x", "y")
     temp<-b[b[,1]==a[i], ]
     for(j in 1:nrow(temp)){
       x<- seq.int(from = (round(temp[j, d])-8), to = (round(temp[j, d])+8), by=1)
       y<- seq.int(from = (round(temp[j, e]-8)), to = (round(temp[j, e])+8), by=1)
       f<- expand.grid(x, y)
       colnames(f)<-c("x", "y")
       mask_df<-rbind(mask_df, f)
       dups<-duplicated(mask_df)
       mask_df<-mask_df[!dups, ]
     }
     mask_list[[i]]<-mask_df[2:nrow(mask_df), ]
  }
  return(mask_list)
}

masks<-mask_gen(images, tiles, 6, 7)

###########################QUANTITY and SIZE ANALYSIS###########################
win_area<-numeric(length = length(images))
feature_density <- numeric(length =length(images))
num_features <- numeric(length =length(images))
feature_area <- numeric(length =length(images))
feature_area_t <- numeric(length =length(images))
rel_feature_area <- numeric(length =length(images))
num_dots<-numeric(length =length(images))
dot_area<-numeric(length =length(images))
dot_area_t<-numeric(length =length(images))
rel_dot_area <-numeric(length = length(images))
dot_density <- numeric(length =length(images))
num_non_dots<-numeric(length =length(images))
non_dot_area<-numeric(length =length(images))
non_dot_area_t<-numeric(length =length(images))
rel_non_dot_area <-numeric(length =length(images))
non_dot_density <-numeric(length =length(images))
num_lipo<-numeric(length =length(images))
lipo_area<-numeric(length =length(images))
lipo_area_t<-numeric(length =length(images))
rel_lipo_area <- numeric(length =length(images))
lipo_density <- numeric(length =length(images))
dx<- character(length = length(images))
for (i in 1:length(images)){
  temp<- data[data$image==images[i], ]
  win_area[i] <- nrow(masks[[i]])
  num_features[i] <- nrow(temp)
  feature_density[i] <- (nrow(temp) / win_area[i])*1000
  feature_area[i] <- mean(temp$area)
  feature_area_t[i] <- sum(temp$area)
  rel_feature_area[i] <- (feature_area_t[i]/win_area[i])*100
  temp1 <- data[data$image==images[i] & data$class=="dots" & 
               data$num_detections>0,]
  dx[i] <- temp1$disease[1]
  num_dots[i] <- nrow(temp1)
  dot_density[i]<- (nrow(temp1)/win_area[i])*1000
  dot_area[i] <- mean(temp1$area)
  dot_area_t[i] <- sum(temp1$area)
  rel_dot_area[i]<-(dot_area_t[i]/win_area[i])*100
  temp2 <- data[data$image==images[i]& data$class=="Non-dots" & 
                data$num_detections>0,]
  num_non_dots[i] <- nrow(temp2)
  non_dot_density[i] <- (nrow(temp2)/win_area[i])*1000
  non_dot_area[i] <- mean(temp2$area)
  non_dot_area_t[i] <- sum(temp2$area)
  rel_non_dot_area[i] <- (non_dot_area_t[i]/win_area[i])*100
  temp3 <- data[data$image==images[i]& data$class=="lipo" & 
                  data$num_detections>0,]
  num_lipo[i] <- nrow(temp3)
  lipo_density[i] <- (nrow(temp3)/win_area[i])*1000
  lipo_area[i]<- mean(temp3$area)
  lipo_area_t[i] <- sum(temp3$area)
  rel_lipo_area[i]<-(lipo_area_t[i]/win_area[i])*100.
}

data_sum<-data.frame(images, dx, win_area, num_features, feature_density, 
                     feature_area, feature_area_t, rel_feature_area, num_dots, 
                     dot_density, dot_area, dot_area_t, rel_dot_area, 
                     num_non_dots, non_dot_density, non_dot_area, non_dot_area_t,
                     rel_non_dot_area, num_lipo, lipo_density, lipo_area, 
                     lipo_area_t, rel_lipo_area)

rm(temp, temp1, temp2, temp3, num_features, feature_density, 
   feature_area, feature_area_t, rel_feature_area, num_dots, 
   dot_density, dot_area, dot_area_t, rel_dot_area, 
   num_non_dots, non_dot_density, non_dot_area, non_dot_area_t,
   rel_non_dot_area, num_lipo, lipo_density, lipo_area, 
   lipo_area_t, rel_lipo_area)

ggplot(data_sum, mapping = aes(x=dx, y=win_area/1000, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Window Area \n (1000 pixels)")+
  xlab("")

wilcox.test(data_sum$win_area~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "win_area"])
median(data_sum[data_sum$dx=="Minimal change", "win_area"])

ggplot(data_sum, mapping = aes(x=dx, y=num_features, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("No. of Features")+
  xlab("")

wilcox.test(data_sum$num_features~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "num_features"])
median(data_sum[data_sum$dx=="Minimal change", "num_features"])

ggplot(data_sum, mapping = aes(x=dx, y=feature_density, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Feature Density \n(per 1000 Pixels)")+
  xlab("")

wilcox.test(data_sum$feature_density~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "feature_density"])
median(data_sum[data_sum$dx=="Minimal change", "feature_density"])

ggplot(data_sum, mapping = aes(x=dx, y=feature_area, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Mean Feature Area \n(pixels)")+
  xlab("")

wilcox.test(data_sum$feature_area~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "feature_area"])
median(data_sum[data_sum$dx=="Minimal change", "feature_area"])

ggplot(data_sum, mapping = aes(x=dx, y=feature_area_t, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Total Feature Area\n(pixels)")+
  xlab("")

wilcox.test(data_sum$feature_area_t~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "feature_area_t"])
median(data_sum[data_sum$dx=="Minimal change", "feature_area_t"])

ggplot(data_sum, mapping = aes(x=dx, y=rel_feature_area, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Relative Feature Area\n(%)")+
  xlab("")+
  ylim(0, 525)

wilcox.test(data_sum$rel_feature_area~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "rel_feature_area"])
median(data_sum[data_sum$dx=="Minimal change", "rel_feature_area"])

ggplot(data_sum, mapping = aes(x=dx, y=num_dots, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("No. of Puncta")+
  xlab("")

wilcox.test(data_sum$num_dots~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "num_dots"])
median(data_sum[data_sum$dx=="Minimal change", "num_dots"])

ggplot(data_sum, mapping = aes(x=dx, y=dot_density, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Puncta Density\n(per 1000 Pixels)")+
  xlab("")

wilcox.test(data_sum$dot_density~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "dot_density"])
median(data_sum[data_sum$dx=="Minimal change", "dot_density"])

ggplot(data_sum, mapping = aes(x=dx, y=dot_area, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Mean Puncta Area\n(pixels)")+
  xlab("")

wilcox.test(data_sum$dot_area~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "dot_area"])
median(data_sum[data_sum$dx=="Minimal change", "dot_area"])

ggplot(data_sum, mapping = aes(x=dx, y=dot_area_t, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Total Puncta Area\n(pixels)")+
  xlab("")

wilcox.test(data_sum$dot_area_t~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "dot_area_t"])
median(data_sum[data_sum$dx=="Minimal change", "dot_area_t"])

ggplot(data_sum, mapping = aes(x=dx, y=rel_dot_area, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Relative Puncta\nArea (%)")+
  xlab("")

wilcox.test(data_sum$rel_dot_area~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "rel_dot_area"])
median(data_sum[data_sum$dx=="Minimal change", "rel_dot_area"])

ggplot(data_sum, mapping = aes(x=dx, y=num_non_dots, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("No. of Non-Puncta\nAreas")+
  xlab("")

wilcox.test(data_sum$num_non_dots~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "num_non_dots"])
median(data_sum[data_sum$dx=="Minimal change", "num_non_dots"])

ggplot(data_sum, mapping = aes(x=dx, y=non_dot_area, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Mean Size Non-\nPuncta Area (pixels)")+
  xlab("")

wilcox.test(data_sum$non_dot_area~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "non_dot_area"])
median(data_sum[data_sum$dx=="Minimal change", "non_dot_area"])

ggplot(data_sum, mapping = aes(x=dx, y=non_dot_area_t/1000, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Total Non-Puncta\nArea (1000 pixels)")+
  xlab("")

wilcox.test(data_sum$non_dot_area_t~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "non_dot_area_t"])
median(data_sum[data_sum$dx=="Minimal change", "non_dot_area_t"])

ggplot(data_sum, mapping = aes(x=dx, y=rel_non_dot_area, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Relative Non-Puncta\nArea (%)")+
  xlab("")+
  ylim(0,121)

wilcox.test(data_sum$rel_non_dot_area~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "rel_non_dot_area"])
median(data_sum[data_sum$dx=="Minimal change", "rel_non_dot_area"])

ggplot(data_sum, mapping = aes(x=dx, y=num_lipo, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("No. of Lipofuscin Dots")+
  xlab("")

wilcox.test(data_sum$num_lipo~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "num_lipo"])
median(data_sum[data_sum$dx=="Minimal change", "num_lipo"])

ggplot(data_sum, mapping = aes(x=dx, y=lipo_area, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Mean Lipofuscin Area\n(pixels)")+
  xlab("")

wilcox.test(data_sum$lipo_area~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "lipo_area"])
median(data_sum[data_sum$dx=="Minimal change", "lipo_area"], na.rm = TRUE)

ggplot(data_sum, mapping = aes(x=dx, y=lipo_area_t, color=dx))+
  geom_boxplot(size=1.8)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  scale_x_discrete(labels=c("MGN", "MCD"))+
  theme(axis.text.x=element_text(size=18, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype="none", color="none")+
  ylab("Total Lipofuscin Area\n(pixels)")+
  xlab("")

wilcox.test(data_sum$lipo_area_t~data_sum$dx)
median(data_sum[data_sum$dx=="Membranous", "lipo_area_t"])
median(data_sum[data_sum$dx=="Minimal change", "lipo_area_t"], na.rm = TRUE)

##################################To Hyperframe#################################
to_hf<-function(a, b, d, e){
  pps<-vector("list", length = length(b))
  dx<-character(length = length(b))
  for (i in 1:length(b)){
    temp<-a[a$image==b[i] & a$parent=="PathAnnotationObject" & a$class==d & 
              a$num_detections>0, ]
    print(nrow(temp))
    print(b[i])
    win<- owin(mask = e[[i]] )
    pps[[i]]<-ppp(temp[,6], temp[,7], window = win)
    dx[i]<-temp[1,14]
  }
  hf<-hyperframe(pps, dx)
  return(hf)
}

dots<-to_hf(data, images, "dots", masks)

##################################Homogeneity###################################
quads<-as.data.frame(matrix(nrow=length(images), ncol =4))
for (i in 1:length(images)){
  quads[i,1]<-dots$V2[i]
  test<-quadrat.test(dots$V1[[i]], method = "MonteCarlo")
  quads[i,2]<-test[[1]]
  quads[i,3]<-test[[3]]
  quads[i,4]<-length(test[[8]])
}
colnames(quads)<-c("Disease", "Chi-Square", "p-Value", "No_Tiles")
write.csv(quads, "MGN_MCD_quadrat_tests.csv")
rm(i, test)

rad_dist<- vector("list", length = length(images))
for (i in 1:nrow(dots)){
  cent<-centroid.owin(Window(dots$V1[[i]]))
  x_dist<-(masks[[i]][,1]-cent[[1]])^2
  y_dist<-(masks[[i]][,2]-cent[[2]])^2
  win_dist<-sqrt(x_dist+y_dist)
  max_dist<- max(win_dist)
  pts<-coords(dots$V1[[i]])
  pts_dist<-sqrt((pts[,1]-cent[[1]])^2+(pts[,2]-cent[[2]])^2)
  pts_rel_dist<- pts_dist/(max_dist)
  rad_dist[[i]]<-data.frame(pts_dist, pts_rel_dist)
}

rm(i, cent, x_dist, y_dist, win_dist, max_dist, pts, pts_dist, pts_rel_dist)

MGN_dist<-rad_dist[dots$V2=="Membranous"]
MGN_pool_dist<-numeric()
for (i in 1:length(MGN_dist)){
  MGN_pool_dist<-c(MGN_pool_dist, MGN_dist[[i]][,2])
}
ggplot(mapping = aes(x=MGN_pool_dist))+
  geom_histogram(fill="black")+
  geom_vline(mapping = aes(xintercept=median(MGN_pool_dist)), color="firebrick3", size=1.5)+
  theme_classic()+
  theme(axis.text.x=element_text(size=15, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 15, colour = "black"),
        axis.title=element_text (size=18), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19),
        plot.title = element_text(size =21, color = "black", vjust = -5, hjust = .05))+
  guides(linetype="none", color="none")+
  ylab("Count")+
  xlab("Relative Centroid Distance")+
  ggtitle("MGN")

MCD_dist<-rad_dist[dots$V2=="Minimal change"]
MCD_pool_dist<-numeric()
for (i in 1:length(MCD_dist)){
  MCD_pool_dist<-c(MCD_pool_dist, MCD_dist[[i]][,2])
}
ggplot(mapping = aes(x=MCD_pool_dist))+
  geom_histogram(fill="black")+
  geom_vline(mapping = aes(xintercept=median(MCD_pool_dist)), color="firebrick3", size=1.5)+
  theme_classic()+
  theme(axis.text.x=element_text(size=15, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 15, colour = "black"),
        axis.title=element_text (size=18), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19),
        plot.title = element_text(size =21, color = "firebrick3", vjust = -5, hjust = .05))+
  guides(linetype="none", color="none")+
  ylab("Count")+
  xlab("Relative Centroid Distance")+
  ggtitle("MCD")

wilcox.test(MGN_pool_dist, MCD_pool_dist)

#########################SPATIAL ANALYSIS IGG MGN vs MCD########################
#Linhom
studpermu.test(dots,  V1~V2, Linhom, correction="translate", use.Tbar = TRUE)

MGN<-dots[dots$V2=="Membranous", ]
MGN_Lest_ls_dots<-lapply(MGN$V1, Linhom, correction="translate")
MGN_Lest_dots_all<-pool.fv(MGN_Lest_ls_dots[[1]], MGN_Lest_ls_dots[[2]], 
                           MGN_Lest_ls_dots[[3]], MGN_Lest_ls_dots[[4]], 
                           MGN_Lest_ls_dots[[5]], MGN_Lest_ls_dots[[6]], 
                           MGN_Lest_ls_dots[[7]], MGN_Lest_ls_dots[[8]])

MCD<-dots[dots$V2 =="Minimal change",]
MCD_Lest_ls_dots<-lapply(MCD$V1, Linhom, correction="translate")
MCD_Lest_dots_all<-pool.fv(MCD_Lest_ls_dots[[1]], MCD_Lest_ls_dots[[2]], 
                           MCD_Lest_ls_dots[[3]], MCD_Lest_ls_dots[[4]], 
                           MCD_Lest_ls_dots[[5]], MCD_Lest_ls_dots[[6]], 
                           MCD_Lest_ls_dots[[7]], MCD_Lest_ls_dots[[8]], 
                           MCD_Lest_ls_dots[[9]], MCD_Lest_ls_dots[[10]], 
                           MCD_Lest_ls_dots[[11]], MCD_Lest_ls_dots[[12]])

r_mgn<-MGN_Lest_dots_all$r
pois_mgn<-MGN_Lest_dots_all$pooltheo
typ<-rep("pois", nrow(MGN_Lest_dots_all))
dx<-rep("MGN", nrow(MGN_Lest_dots_all))
pois_mgn<-cbind(r_mgn, pois_mgn, typ, dx)
colnames(pois_mgn)<-c("r", "value", "type", "dx")
est_mgn<-MGN_Lest_dots_all$pooltrans
typ<-rep("est", nrow(MGN_Lest_dots_all))
dx<-rep("MGN", nrow(MGN_Lest_dots_all))
est_mgn<-cbind(r_mgn, est_mgn, typ, dx)
colnames(est_mgn)<-c("r", "value", "type", "dx")

MGN_plot<-as.data.frame(rbind(pois_mgn, est_mgn))
MGN_plot[,1]<- as.numeric(MGN_plot[,1])
MGN_plot[,2]<- as.numeric(MGN_plot[,2])
MGN_plot[,3]<- as.factor(MGN_plot[,3])
MGN_plot[,4]<- as.factor(MGN_plot[,4])

ggplot(MGN_plot, mapping = aes(x=r, y=value, linetype=type))+
  geom_line(size=1.4)+
  geom_ribbon(aes(ymin=c(MGN_Lest_dots_all$lotrans,MGN_Lest_dots_all$lotrans), 
                  ymax=c(MGN_Lest_dots_all$hitrans, MGN_Lest_dots_all$hitrans)), 
              alpha=0.2)+
  theme_classic()+
  theme(axis.text.x=element_text(size=18, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  scale_color_manual(values = "black")+
  theme(axis.text.x=element_text(size=18), axis.text.y = element_text(size = 18),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype=FALSE)+
  ylab("L(r)")

r_mcd<-MCD_Lest_dots_all$r
pois_mcd<-MCD_Lest_dots_all$pooltheo
typ<-rep("pois", nrow(MCD_Lest_dots_all))
dx<-rep("MCD", nrow(MCD_Lest_dots_all))
pois_mcd<-cbind(r_mcd, pois_mcd, typ, dx)
colnames(pois_mcd)<-c("r", "value", "type", "dx")
est_mcd<-MCD_Lest_dots_all$pooltrans
typ<-rep("est", nrow(MCD_Lest_dots_all))
dx<-rep("MCD", nrow(MCD_Lest_dots_all))
est_mcd<-cbind(r_mcd, est_mcd, typ, dx)
colnames(est_mcd)<-c("r", "value", "type", "dx")

MCD_plot<-as.data.frame(rbind(pois_mcd, est_mcd))
MCD_plot[,1]<- as.numeric(MCD_plot[,1])
MCD_plot[,2]<- as.numeric(MCD_plot[,2])
MCD_plot[,3]<- as.factor(MCD_plot[,3])
MCD_plot[,4]<- as.factor(MCD_plot[,4])

ggplot(MCD_plot, mapping = aes(x=r, y=value, linetype=type, color=type))+
  geom_line(size=1.4)+
  geom_ribbon(aes(ymin=c(MCD_Lest_dots_all$lotrans,MCD_Lest_dots_all$lotrans), 
                  ymax=c(MCD_Lest_dots_all$hitrans, MCD_Lest_dots_all$hitrans)), 
              alpha=0.2)+
  theme_classic()+
  theme(axis.text.x=element_text(size=18, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  scale_color_manual(values = c("firebrick3", "black"))+
  theme(axis.text.x=element_text(size=18), axis.text.y = element_text(size = 18),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype=FALSE, color=FALSE )+
  ylab("L(r)")

comp_plot<-rbind(MGN_plot, MCD_plot)

ggplot(comp_plot, mapping = aes(x=r, y=value, color=dx, linetype=type))+
  geom_line(size=1.5)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  theme(axis.text.x=element_text(size=18, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype=FALSE)+
  labs(color= "Disease")+
  ylab("L(r)")

#PCF
studpermu.test(dots,  V1~V2, pcf, method="d", divisor="d", correction="translate", use.Tbar = TRUE)

MGN<-dots[dots$V2=="Membranous", ]
MGN_Lest_ls_dots<-lapply(MGN$V1, pcf, method="d", divisor="d", correction="translate")
MGN_Lest_dots_all<-pool.fv(MGN_Lest_ls_dots[[1]], MGN_Lest_ls_dots[[2]], 
                           MGN_Lest_ls_dots[[3]], MGN_Lest_ls_dots[[4]], 
                           MGN_Lest_ls_dots[[5]], MGN_Lest_ls_dots[[6]], 
                           MGN_Lest_ls_dots[[7]], MGN_Lest_ls_dots[[8]])

MCD<-dots[dots$V2 =="Minimal change",]
MCD_Lest_ls_dots<-lapply(MCD$V1, pcf, method="d", divisor="d", correction="translate")
MCD_Lest_dots_all<-pool.fv(MCD_Lest_ls_dots[[1]], MCD_Lest_ls_dots[[2]], 
                           MCD_Lest_ls_dots[[3]], MCD_Lest_ls_dots[[4]], 
                           MCD_Lest_ls_dots[[5]], MCD_Lest_ls_dots[[6]], 
                           MCD_Lest_ls_dots[[7]], MCD_Lest_ls_dots[[8]], 
                           MCD_Lest_ls_dots[[9]], MCD_Lest_ls_dots[[10]], 
                           MCD_Lest_ls_dots[[11]], MCD_Lest_ls_dots[[12]])

r_mgn<-MGN_Lest_dots_all$r
pois_mgn<-MGN_Lest_dots_all$pooltheo
typ<-rep("pois", nrow(MGN_Lest_dots_all))
dx<-rep("MGN", nrow(MGN_Lest_dots_all))
pois_mgn<-cbind(r_mgn, pois_mgn, typ, dx)
colnames(pois_mgn)<-c("r", "value", "type", "dx")
est_mgn<-MGN_Lest_dots_all$pooltrans
typ<-rep("est", nrow(MGN_Lest_dots_all))
dx<-rep("MGN", nrow(MGN_Lest_dots_all))
est_mgn<-cbind(r_mgn, est_mgn, typ, dx)
colnames(est_mgn)<-c("r", "value", "type", "dx")

MGN_plot<-as.data.frame(rbind(pois_mgn, est_mgn))
MGN_plot[,1]<- as.numeric(MGN_plot[,1])
MGN_plot[,2]<- as.numeric(MGN_plot[,2])
MGN_plot[,3]<- as.factor(MGN_plot[,3])
MGN_plot[,4]<- as.factor(MGN_plot[,4])

r_mcd<-MCD_Lest_dots_all$r
pois_mcd<-MCD_Lest_dots_all$pooltheo
typ<-rep("pois", nrow(MCD_Lest_dots_all))
dx<-rep("MCD", nrow(MCD_Lest_dots_all))
pois_mcd<-cbind(r_mcd, pois_mcd, typ, dx)
colnames(pois_mcd)<-c("r", "value", "type", "dx")
est_mcd<-MCD_Lest_dots_all$pooltrans
typ<-rep("est", nrow(MCD_Lest_dots_all))
dx<-rep("MCD", nrow(MCD_Lest_dots_all))
est_mcd<-cbind(r_mcd, est_mcd, typ, dx)
colnames(est_mcd)<-c("r", "value", "type", "dx")

MCD_plot<-as.data.frame(rbind(pois_mcd, est_mcd))
MCD_plot[,1]<- as.numeric(MCD_plot[,1])
MCD_plot[,2]<- as.numeric(MCD_plot[,2])
MCD_plot[,3]<- as.factor(MCD_plot[,3])
MCD_plot[,4]<- as.factor(MCD_plot[,4])

comp_plot<-rbind(MGN_plot, MCD_plot)
rm(typ, dx, r_mcd, r_mgn, pois_mcd, pois_mgn, est_mcd, est_mgn, MGN_plot, MCD_plot)

ggplot(comp_plot, mapping = aes(x=r, y=value, color=dx, linetype=type))+
  geom_line(size=1.5)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  theme(axis.text.x=element_text(size=18, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype=FALSE)+
  labs(color= "Disease")+
  ylab("Pair Correlation (r)")

#G inhom Nearest Neighbor Function
studpermu.test(dots,  V1~V2, Ginhom, correction="translate", use.Tbar = TRUE)

MGN<-dots[dots$V2=="Membranous", ]
MGN_Lest_ls_dots<-lapply(MGN$V1, Ginhom, correction="translate")
MGN_Lest_dots_all<-pool.fv(MGN_Lest_ls_dots[[1]], MGN_Lest_ls_dots[[2]], 
                           MGN_Lest_ls_dots[[3]], MGN_Lest_ls_dots[[4]], 
                           MGN_Lest_ls_dots[[5]], MGN_Lest_ls_dots[[6]], 
                           MGN_Lest_ls_dots[[7]], MGN_Lest_ls_dots[[8]])

MCD<-dots[dots$V2 =="Minimal change",]
MCD_Lest_ls_dots<-lapply(MCD$V1, Ginhom, correction="translate")
MCD_Lest_dots_all<-pool.fv(MCD_Lest_ls_dots[[1]], MCD_Lest_ls_dots[[2]], 
                           MCD_Lest_ls_dots[[3]], MCD_Lest_ls_dots[[4]], 
                           MCD_Lest_ls_dots[[5]], MCD_Lest_ls_dots[[6]], 
                           MCD_Lest_ls_dots[[7]], MCD_Lest_ls_dots[[8]], 
                           MCD_Lest_ls_dots[[9]], MCD_Lest_ls_dots[[10]], 
                           MCD_Lest_ls_dots[[11]], MCD_Lest_ls_dots[[12]])

r_mgn<-MGN_Lest_dots_all$r
pois_mgn<-MGN_Lest_dots_all$pooltheo
typ<-rep("pois", nrow(MGN_Lest_dots_all))
dx<-rep("MGN", nrow(MGN_Lest_dots_all))
pois_mgn<-cbind(r_mgn, pois_mgn, typ, dx)
colnames(pois_mgn)<-c("r", "value", "type", "dx")
est_mgn<-MGN_Lest_dots_all$poolbord
typ<-rep("est", nrow(MGN_Lest_dots_all))
dx<-rep("MGN", nrow(MGN_Lest_dots_all))
est_mgn<-cbind(r_mgn, est_mgn, typ, dx)
colnames(est_mgn)<-c("r", "value", "type", "dx")

MGN_plot<-as.data.frame(rbind(pois_mgn, est_mgn))
MGN_plot[,1]<- as.numeric(MGN_plot[,1])
MGN_plot[,2]<- as.numeric(MGN_plot[,2])
MGN_plot[,3]<- as.factor(MGN_plot[,3])
MGN_plot[,4]<- as.factor(MGN_plot[,4])

r_mcd<-MCD_Lest_dots_all$r
pois_mcd<-MCD_Lest_dots_all$pooltheo
typ<-rep("pois", nrow(MCD_Lest_dots_all))
dx<-rep("MCD", nrow(MCD_Lest_dots_all))
pois_mcd<-cbind(r_mcd, pois_mcd, typ, dx)
colnames(pois_mcd)<-c("r", "value", "type", "dx")
est_mcd<-MCD_Lest_dots_all$poolbord
typ<-rep("est", nrow(MCD_Lest_dots_all))
dx<-rep("MCD", nrow(MCD_Lest_dots_all))
est_mcd<-cbind(r_mcd, est_mcd, typ, dx)
colnames(est_mcd)<-c("r", "value", "type", "dx")

MCD_plot<-as.data.frame(rbind(pois_mcd, est_mcd))
MCD_plot[,1]<- as.numeric(MCD_plot[,1])
MCD_plot[,2]<- as.numeric(MCD_plot[,2])
MCD_plot[,3]<- as.factor(MCD_plot[,3])
MCD_plot[,4]<- as.factor(MCD_plot[,4])

comp_plot<-rbind(MGN_plot, MCD_plot)
rm(typ, dx, r_mcd, r_mgn, pois_mcd, pois_mgn, est_mcd, est_mgn, MGN_plot, MCD_plot)

ggplot(comp_plot, mapping = aes(x=r, y=value, color=dx, linetype=type))+
  geom_line(size=1.5)+
  theme_classic()+
  scale_color_manual(values = c("black", "firebrick3"))+
  theme(axis.text.x=element_text(size=18, colour = "black"), 
        axis.text.y = element_text(size = 18, colour = "black"),
        axis.title=element_text (size=21), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19))+
  guides(linetype=FALSE)+
  labs(color= "Disease")+
  ylab("G(r)")

############################## QA for AI of Dot space ##########################

XXXXXXX_val<-read.table(
  "FILE_PATH",
  header = TRUE, stringsAsFactors = FALSE, check.names = FALSE, sep = "\t")
XXXXXXX_val<-read.table(
  "FILE_PATH",
  header = TRUE, stringsAsFactors = FALSE, check.names = FALSE, sep = "\t")
XXXXXXXX_IgG_val<-read.table(
  "FILE_PATH",
  header = TRUE, stringsAsFactors = FALSE, check.names = FALSE, sep = "\t")
XXXXXXXX_val<-read.table(
  "FILE_PATH",
  header = TRUE, stringsAsFactors = FALSE, check.names = FALSE, sep = "\t")
XXXXXXXX_IgG_val<-read.table(
  "FILE_PATH",
  header = TRUE, stringsAsFactors = FALSE, check.names = FALSE, sep = "\t")
XXXXXXX_val<-read.table(
  "FILE_PATH",
  header = TRUE, stringsAsFactors = FALSE, check.names = FALSE, sep = "\t")

XXXXX_val_win1<- disc(radius = 75.2396, centre = c(1320.5, 464))
XXXXX_val_win2<- disc(radius = 96.7594, centre = c(921.5, 638))

XXXXX_val_win1<-disc(radius = 63.4521, centre = c(1398.5, 981))
XXXXX_val_win2<-disc(radius = 100.2348, centre = c(1487.5, 505))

XXXXX_val_win1<- disc(radius = 86.5506, centre = c(955.5, 702.5))
XXXXX_val_win2<- disc(radius = 74.8031, centre = c(764.5, 416.5))

XXXXX_val_win1<- disc(radius = 53.2068, centre = c(811.5, 283.5))
XXXXX_val_win2<- disc(radius = 56.5147, centre = c(1103.5, 658))

XXXXX_val_win1<- disc(radius = 39.8925, centre = c(985, 487))
XXXXX_val_win2<- disc(radius = 42.0940, centre = c(997, 730.5))

XXXXX_IgG_val1<-ppp(XXXXX_IgG_val$x, XXXXXX_IgG_val$y, 
                       XXXXXX_val_win1)
XXXXXX_IgG_val2<-ppp(XXXXXX_IgG_val$x, XXXXXXX_IgG_val$y, 
                       XXXXXXX_val_win2)
XXXXXX_IgG_test1<-dots[4,1, drop=TRUE]
Window(XXXXXXXX_IgG_test1)<-XXXXXX_val_win1
XXXXXX_IgG_test2<-dots[4,1, drop=TRUE]
Window(XXXXXXX_IgG_test2)<-XXXXXXXX_val_win2

XXXXXX_IgG_val1<-ppp(XXXXXXX_IgG_val$x, XXXXXXX_IgG_val$y, 
                        XXXXXXXX_val_win1)
XXXXXX_IgG_val2<-ppp(XXXXXX_IgG_val$x, XXXXXX_IgG_val$y, 
                        XXXXXX_val_win2)
XXXXXXX_IgG_test1<-dots[13, 1, drop=TRUE]
Window(XXXXX_IgG_test1)<-XXXXXX_val_win1
XXXXXXXX_IgG_test2<-dots[13, 1, drop=TRUE]
Window(XXXXXXX_IgG_test2)<-XXXXXXX_val_win2


XXXXXX_IgG_val1<-ppp(XXXXXX_IgG_val$x, XXXXXX_IgG_val$y, 
                        XXXXXX_val_win1)
XXXXXX_IgG_val2<-ppp(XXXXXX_IgG_val$x, XXXXXXX_IgG_val$y, 
                        XXXXXXX_val_win2)
XXXXXX_IgG_test1<-dots[12, 1, drop=TRUE]
Window(XXXXXXX_IgG_test1)<-XXXXXXX_val_win1
XXXXXX_IgG_test2<-dots[12, 1, drop=TRUE]
Window(XXXXXX_IgG_test2)<-XXXXXXXX_val_win2
XXXXXX_IgG_val1<-ppp(XXXXXXX_IgG_val$x, XXXXXXXX_IgG_val$y, 
                         XXXXXXX_val_win1)
XXXXXXXX_IgG_val2<-ppp(XXXXXXX_IgG_val$x, XXXXXXX_IgG_val$y, 
                         XXXXXXX_val_win2)
XXXXXX_IgG_test1<-dots[15, 1, drop=TRUE]
Window(XXXXXX_IgG_test1)<-XXXXXX_val_win1
XXXXXXXX_IgG_test2<-dots[15, 1, drop=TRUE]
Window(XXXXXXXX_IgG_test2)<-XXXXXXX_val_win2

XXXXXX_IgG_val1<-ppp(XXXXXX_IgG_val$x, XXXXX_IgG_val$y, 
                         XXXXXX_val_win1)
XXXXXX_IgG_val2<-ppp(XXXXXXXX_IgG_val$x, XXXXXXX_IgG_val$y, 
                         XXXXXXXX_val_win2)
XXXX_IgG_test1<-dots[16, 1, drop=TRUE]
Window(XXXXXX_IgG_test1)<-XXXXXX_val_win1
XXXXX_IgG_test2<-dots[16, 1, drop=TRUE]
Window(XXXXXXX_IgG_test2)<-XXXXXX_val_win2

val_comp<-function(test, val){
  test_coord<-coords(test)
  val_coord<-coords(val)
  dist_mat<- as.data.frame(matrix(nrow = (nrow(test_coord)*nrow(val_coord)), ncol=3))
  k=1
  for (i in 1:nrow(test_coord)){
    test_point<-test_coord[i,]
    for (j in 1:nrow(val_coord)){
      dist_mat[k,1]<-i
      dist_mat[k,2]<-j
      dist_mat[k,3]<-sqrt((test_point[1]-val_coord[j,1])^2+(test_point[2]-val_coord[j,2])^2)
      k=k+1
    }
  }
  nn_dist<-as.data.frame(matrix(nrow = nrow(val_coord), ncol = 3))
  for (i in 1:nrow(val_coord)){
    temp<-dist_mat[dist_mat[,2]==i, ]
    nn_dist[i,]<-temp[temp[,3]==min(temp[,3]), ]
  }
  dist_mat<-dist_mat[dist_mat[,3]<7,]
  matched_points<-as.data.frame(matrix(ncol = 3))
  temp<-dist_mat
  while (nrow(temp)>0){
    min_dist<- temp[temp[,3]==min(temp[,3]),]
    matched_points<-rbind(matched_points, min_dist)
    temp<-temp[temp[,1]!=min_dist[1,1],]
    temp<-temp[temp[,2]!=min_dist[1,2],]
  }
  master<-numeric(length = 6)
  master[1]<-nrow(test_coord)
  master[2]<-nrow(val_coord)
  master[3]<-nrow(matched_points)-1
  master[4]<-nrow(dist_mat)
  master[5]<-master[1]-master[3]
  master[6]<-master[2]-master[3]
  names(master)<-c("test_pts", "val_pts", "matched_pts", "matched_area", "unmatch_pts_test", "unmatched_pts_val")
  results<-c(list(nn_dist), list(master), list(dist_mat))
  return(results)
}

val_comp_XXXX_1<-val_comp(test = XXXX_IgG_test1, val = XXXX_IgG_val1)
val_comp_XXXX_2<-val_comp(test = XXXX_IgG_test2, val = XXXX_IgG_val2)

val_comp_XXXX_1<-val_comp(test = XXXX_IgG_test1, val = XXXX_IgG_val1)
val_comp_XXXX_2<-val_comp(test = XXXX_IgG_test2, val = XXXX_IgG_val2)

val_comp_XXXX_1<-val_comp(test = XXXX_IgG_test1, val = XXXX_IgG_val1)
val_comp_XXXX_2<-val_comp(test = XXXX_IgG_test2, val = XXXX_IgG_val2)

val_comp_XXXX_1<-val_comp(test = XXXX_IgG_test1, val = XXXX_IgG_val1)
val_comp_XXXX_2<-val_comp(test = XXXX_IgG_test2, val = XXXX_IgG_val2)

val_comp_XXXX_1<-val_comp(test = XXXX_IgG_test1, val = XXXX_IgG_val1)
val_comp_XXXX_2<-val_comp(test = XXXX_IgG_test2, val = XXXX_IgG_val2)

combined_nn_dist<- rbind(val_comp_XXXX_1[[1]], val_comp_XXXX_2[[1]], 
                         val_comp_XXXX_1[[1]], val_comp_XXXXX_2[[1]], 
                         val_comp_XXXX_1[[1]], val_comp_XXXX_2[[1]], 
                         val_comp_XXXX_1[[1]], val_comp_XXXXX_2[[1]], 
                         val_comp_XXXX_1[[1]], val_comp_XXXX_2[[1]])
ggplot(combined_nn_dist, mapping = aes(x=combined_nn_dist$V3))+
  geom_histogram(fill="black")+
  theme_classic()+
  theme(axis.text.x=element_text(size=15, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 15, colour = "black"),
        axis.title=element_text (size=18), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19),
        plot.title = element_text(size =21, color = "black", vjust = -5, hjust = .05))+
  guides(linetype="none", color="none")+
  ylab("Count")+
  xlab("Distance (Pixels)")

combined_dist_mat<-rbind(val_comp_XXXX_1[[3]], val_comp_XXXX_2[[3]], 
                         val_comp_XXXX_1[[3]], val_comp_XXXX_2[[3]], 
                         val_comp_XXXX_1[[3]], val_comp_XXXX_2[[3]], 
                         val_comp_XXXX_1[[3]], val_comp_XXXX_2[[3]], 
                         val_comp_XXXX_1[[3]], val_comp_XXXX_2[[3]])
ggplot(combined_dist_mat, mapping = aes(x=combined_dist_mat$V3))+
  geom_histogram(fill="black")+
  theme_classic()+
  theme(axis.text.x=element_text(size=15, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 15, colour = "black"),
        axis.title=element_text (size=18), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19),
        plot.title = element_text(size =21, color = "black", vjust = -5, hjust = .05))+
  guides(linetype="none", color="none")+
  ylab("Count")+
  xlab("Distance (Pixels)")

combined_val_comp<-t(data.frame(val_comp_XXXX_1[[2]], val_comp_XXXX_2[[2]], 
                              val_comp_XXXX_1[[2]], val_comp_XXXX_2[[2]], 
                              val_comp_XXXX_1[[2]], val_comp_XXXX_2[[2]], 
                              val_comp_XXXX_1[[2]], val_comp_XXXX_2[[2]], 
                              val_comp_XXXX_1[[2]], val_comp_XXXX_2[[2]]))

sensitivity_OA<-numeric(length = 3)
sensitivity_OA[1]<-sum(combined_val_comp[1:6,4])/sum(combined_val_comp[1:6,2])
sensitivity_OA[2]<-sum(combined_val_comp[7:10,4])/sum(combined_val_comp[7:10,2])
sensitivity_OA[3]<-sum(combined_val_comp[,4])/sum(combined_val_comp[,2])
dx<-c("MCD", "MGN", "Overall")
sensitivity_OA<-as.data.frame(cbind(dx, sensitivity_OA))
sensitivity_OA[,2]<-as.numeric(sensitivity_OA[,2])
ggplot(mapping = aes(x=sensitivity_OA$dx, y=sensitivity_OA$sensitivity_OA))+
  geom_col(fill="black")+
  theme_classic()+
  theme(axis.text.x=element_text(size=15, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 15, colour = "black"),
        axis.title=element_text (size=18), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19),
        plot.title = element_text(size =21, color = "black", vjust = -5, hjust = .05))+
  guides(linetype="none", color="none")+
  ylab("Sensitivity")+
  xlab("Dx")

sensitivity_sw<-numeric(length = 10)
for (i in 1:10){
  sensitivity_sw[i]<-combined_val_comp[i,4]/combined_val_comp[i,2]
}

wilcox.test(sensitivity_sw[1:6], sensitivity_sw[7:10]) #p=0.3524

FP_OA<-numeric(length = 3)
FP_OA[1]<-sum(combined_val_comp[1:6,5])/sum(combined_val_comp[1:6,1])
FP_OA[2]<-sum(combined_val_comp[7:10,5])/sum(combined_val_comp[7:10,1])
FP_OA[3]<-sum(combined_val_comp[,5])/sum(combined_val_comp[,1])
dx<-c("MCD", "MGN", "Overall")
FP_OA<-as.data.frame(cbind(dx, FP_OA))
FP_OA[,2]<-as.numeric(FP_OA[,2])
ggplot(mapping = aes(x=FP_OA$dx, y=FP_OA$FP_OA))+
  geom_col(fill="black")+
  theme_classic()+
  theme(axis.text.x=element_text(size=15, angle = 45, hjust = 1.0, colour = "black"), 
        axis.text.y = element_text(size = 15, colour = "black"),
        axis.title=element_text (size=18), legend.text=element_text(size=17), 
        legend.title = element_text(size = 19),
        plot.title = element_text(size =21, color = "black", vjust = -5, hjust = .05))+
  guides(linetype="none", color="none")+
  ylab("False Positive Rate")+
  xlab("Dx")

FP_sw<-numeric(length = 10)
for (i in 1:10){
  FP_sw[i]<-combined_val_comp[i,5]/combined_val_comp[i,1]
}

wilcox.test(FP_sw[1:6], FP_sw[7:10]) #p=0.161

for(i in 1:length(test)){
  temp<-which(data$image==test[i])
  de_id[temp]<-test2[i]
  temp2<-which(tiles$Image==test[i])
  de_id2[temp2]<-test2[i]
}
data2<-cbind(data[,1], de_id, data[,2:ncol(data)])
write.csv(data2, "raw_data_supp.csv")
tiles2<-cbind(tiles[,1], de_id2, tiles[2:ncol(tiles)])
write.csv(tiles2, "tiles_supp.csv")